/*
 * pt.cpp
 *
 *  Created on: Nov 20, 2010
 *      Author: snirgaz
 */

#include "ParallelTempring.h"
using namespace std;
ParallelTempring::ParallelTempring(QmcParams &SP, PTParams &PP) :
		QmcParams(SP), PTParams(PP) {
	DoubleVec::iterator ir;
	IntVec::iterator ip;
	StringVec::iterator id, ii;
	int seed = QmcParams::seed, i;
	QmcParams CP = static_cast<QmcParams>(*this);
	confData cd;
	int max_delta = 3;
	histUp_.resize(PTParams::numOfConf, 0);
	histDown_.resize(PTParams::numOfConf, 0);
	ptRate_.resize(2, 0);
	std::fill(ptRate_.begin(), ptRate_.end(), 0);
	this->initializeConfHistogram();
	OutPutfilePT_.open((PTParams::PTOutFile).c_str(), std::ios::out);
	// LAd rs from a file
	if (PTParams::initType == "FILE_CONTINUE_RUN") {
		// Load stream position
		String file_name = PTParams::initPath + "/randomStreamPT";
		ptStream_.loadStream(file_name);
		CP.init = true;
		CP.initOnlyWLData = false;
	}
	if (PTParams::initType == "FILE_FIRST_RUN") {
		CP.init = false;
		CP.initOnlyWLData = true;
	}
	if ((PTParams::initType == "FILE_CONTINUE_RUN")
			|| (PTParams::initType == "FILE_FIRST_RUN")) {
		String file_name = PTParams::initPath + "RSS_Final.h5";
		String dsName = "/DS1";
		hdf5_rw::read_double(file_name, dsName, &(PTParams::RSS[0]));
	} else {
		ptStream_.initStream(BRNG + seed++);
		CP.init = false;
		CP.initOnlyWLData = false;
	}
	for (ir = PTParams::RSS.begin(), id = PTParams::outDirs.begin(), ii =
			PTParams::initDirs.begin(), i = 0; ir != PTParams::RSS.end();
			ir++, id++, i++) {
		// Different Seeds
		CP.seed = seed++;
#ifdef ROT
		CP.lambda = 1
		/ (PI * double(QmcParams::NumP) * pow(*ir, 2) / (sqrt(3) / 2.));
#else
		CP.lambda = 1 / (PI * double(QmcParams::NumP) * pow(*ir, 2));
#endif
		CP.epsilon = (1 / ((double) (CP.M))) * CP.beta;
		CP.sigma = sqrt(2 * CP.lambda * CP.epsilon);
		Sigmas_.push_back(0.5 * pow(CP.sigma, -2));
		Force_Factor.push_back(
				(1. / 9.) * pow(CP.epsilon, 3) * (2. * CP.lambda));
		CP.outPath = *id;
		if (CP.init || CP.initOnlyWLData)
			CP.initPath = (*ii++);
		CP.confNum = i;
		confs_.push_back(new Qmc(CP));
	}
	for (i = 0; i < PTParams::numOfConf; i++) {
		cd.confNum = i;
		cd.confDir = UP;
		confsData_.push_back(cd);
	}
	// PT random Number Initialize
	jump_.init(QmcParams::rndBufferSize, &ptStream_, 1, PTParams::maxDelta + 1);
	upDown_.init(QmcParams::rndBufferSize, &ptStream_, 0, 2);
	conf_.init(QmcParams::rndBufferSize, &ptStream_, 0, PTParams::numOfConf);
	uniZeroOne_.init(QmcParams::rndBufferSize, &ptStream_, 0, 1);
	// Timing
	ippGetCpuFreqMhz(&mhz_);
	stopTimeEnd_ = false;
}
ParallelTempring::~ParallelTempring() {
}
void ParallelTempring::thermalize() {
#pragma omp parallel for schedule(static)
	for (int i = 0; i < PTParams::numOfConf; i++) {
		confs_[i].thermalize();
	}
}
void ParallelTempring::initializeFromFile() {
	for (int i = 0; i < PTParams::numOfConf; i++)
		confs_[i].initializeFromFile();
}
void ParallelTempring::initializeWLDataFromFile() {
	for (int i = 0; i < PTParams::numOfConf; i++)
		confs_[i].initialize();
}

void ParallelTempring::initializeSolid() {
	for (int i = 0; i < PTParams::numOfConf; i++)
		(confs_[i]).initialize();
}

void ParallelTempring::run() {
	stopTimeEnd_ = false;
	bool stopMeasure = false;
	long int run_num = 0;
	long int measureNum = 0;
	Ipp64u start, end;
	if (PTParams::initType == "FILE_CONTINUE_RUN") {
		OutPutfilePT_ << "Continue Last Run " << std::endl;
		this->initializeFromFile();
		this->maxEnergyDiff();
	}
	if (PTParams::initType == "FILE_FIRST_RUN") {
		OutPutfilePT_ << "First Run - Initialize From File " << std::endl;
		this->initializeWLDataFromFile();
#ifdef PT_PAR
		this->adjustLamabdas();
#endif
		this->maxEnergyDiff();
	}
	if (PTParams::initType == "SOLID_FIRST_RUN") {
		OutPutfilePT_ << "First Run - Initialize As Solid " << std::endl;
		this->initializeSolid();
#ifdef PT_PAR
		this->adjustLamabdas();
#endif
		this->maxEnergyDiff();
	}
	int maxSweepRate = this->getMaxRate();
	gettimeofday(&startAll_, NULL);
	std::for_each(confs_.begin(), confs_.end(),
			[](Qmc &conf) {conf.initializeCountOutput();});
	while (!stopMeasure) {
		// run
		//start = ippGetCpuClocks();
		measureNum++;
		for (int j = 0; j < PTParams::numSplit; j++) {
			this->runN(maxSweepRate);
		}
#ifdef PT_PAR
//		OutPutfilePT_ << "Energy Diff after run";
//		this->maxEnergyDiff();
		for (int i = 0; i < PTParams::numPar; i++)
		this->swapConfMove();
//		OutPutfilePT_ << "Energy Diff after swap";
//		this->maxEnergyDiff();
//
//	}
#endif
//		this->swapConfMove();
		if (!((run_num++) % PTParams::numRunsPrint)) {
			this->printStatus();
			this->updateConfHistogram();
		}
		this->ptMeasure();
#ifdef SPLIT_SIM
		stopMeasure = ((this->ptMeasureDone()) || stopTimeEnd_);
#else
		stopMeasure = this->ptMeasureDone();
#endif
	}
	this->maxEnergyDiff();
}
void ParallelTempring::runN(int n) {
#pragma omp parallel for schedule(static)
	for (int i = 0; i < PTParams::numOfConf; i++) {
		confs_[i].runN(n);
	}
}
void ParallelTempring::ptMeasure() {
#pragma omp parallel for
	for (int i = 0; i < PTParams::numOfConf; i++) {
		confs_[i].performMeasure();
	}
}
int ParallelTempring::getMaxRate() {
	IntVec rates;
	for (int i = 0; i < PTParams::numOfConf; i++) {
		rates.push_back(confs_[i].getSweepRate());
	}
	return *(std::max_element(rates.begin(), rates.end()));
}
void ParallelTempring::swapConfMove() {
	int conf1, conf2, up_down;
	double dS, p;
	int delta, jump;
// Draw conf1 , conf2
	conf1 = conf_.getSingle();
	if (conf1 == 0)
		up_down = 1;
	else if (conf1 == PTParams::numOfConf - 1)
		up_down = 0;
	else {
		up_down = upDown_.getSingle();
	}
	jump = jump_.getSingle();
	if (up_down) {
		delta = fmin(jump, PTParams::numOfConf - 1 - conf1);
		conf2 = conf1 + delta;
	} else {
		delta = fmin(jump, conf1);
		conf2 = conf1 - delta;
	}
	dS = this->calcDS(conf1, conf2);
	p = exp(dS);
	if (this->accept(p)) {
		this->swapConf(conf1, conf2);
		ptRate_[1] += 1;
	} else
		ptRate_[0] += 1;
}
double ParallelTempring::calcDS(int conf1, int conf2) {
	double l1, l2, f1, f2, Ek1, Ek2, dS1, dS2, ratio;
	WLEnergyData ED1, ED2;
	l1 = Sigmas_[conf1];
	l2 = Sigmas_[conf2];
	f1 = Force_Factor[conf1];
	f2 = Force_Factor[conf2];
	Ek1 = (confs_[conf1]).getEK();
	Ek2 = (confs_[conf2]).getEK();
	ED1 = (confs_[conf1]).getED();
	ED2 = (confs_[conf2]).getED();
	assert((ED1.energyF_ > 0) && (ED2.energyF_ > 0));
	dS1 = (l1 - l2) * (Ek1 - Ek2);
	dS2 = (f1 - f2) * (ED1.getEnergyF() - ED2.getEnergyF());
	return dS1 + dS2;
}
bool ParallelTempring::accept(double p) {
	double draw;
	if (p < 1) {
		// Draw uniform [0,1]
		draw = uniZeroOne_.getSingle();
		// Check if below boltzman weight
		if (draw > p)
			return false;
	}
	return true;
}
void ParallelTempring::swapConf(int conf1, int conf2) {
	Qmc &qmcConf1 = confs_[conf1], &qmcConf2 = confs_[conf2];
	qmcConf1.swapWL(qmcConf2);
	this->updateConfVec(conf1, conf2);
}
bool ParallelTempring::ptMeasureDone() {
	stopMeasureEnd_ = true;
	for (int i = 0; i < PTParams::numOfConf; i++) {
		stopMeasureEnd_ = stopMeasureEnd_ && (confs_[i].measureDone());
	}
	if (stopMeasureEnd_)
		OutPutfilePT_ << "\n Measure Ended" << std::flush;
	return stopMeasureEnd_;
}
void ParallelTempring::dump() {
	std::string param, file_name;
	for (int i = 0; i < PTParams::numOfConf; i++)
		confs_[i].dump();
	file_name = PTParams::outPathPT;
// Write Parameters to XML file
	file_name.append("PT_Output.xml");
	xmlWriter xml_output(file_name);
	xml_output.newElement(param = "PT_PARAMETERS");
	xml_output.setDataElement(param = "ACCEPT_PT", (double) ptRate_[1]);
	xml_output.setDataElement(param = "REJECT_PT", (double) ptRate_[0]);
	xml_output.setDataElement(param = "MEASURE_DONE", stopMeasureEnd_);
	xml_output.setDataElement(param = "RUN_NUM", PTParams::runNum);
	xml_output.writeToFile();
	this->dumpHistograms(-1);
// Save stream position
	file_name = PTParams::outPathPT + "/randomStreamPT";
	ptStream_.saveStream(file_name);
}
void ParallelTempring::dumpConf() {
	for (int i = 0; i < PTParams::numOfConf; i++) {
		confs_[i].dump();
	}
}
void ParallelTempring::printConf() {
	for (int i = 0; i < PTParams::numOfConf; i++)
		OutPutfilePT_ << confsData_[i].confNum << "-";
}
int ParallelTempring::countClosed() {
	int sum = 0;
	for (int i = 0; i < PTParams::numOfConf; i++)
		if (confs_[i].isClosed())
			sum++;
	return sum;
}
void ParallelTempring::updateConfVec(int conf1, int conf2) {
	confData pos1, pos2;
	pos1 = confsData_[conf1];
	pos2 = confsData_[conf2];
	confsData_[conf1] = pos2;
	confsData_[conf2] = pos1;
	if ((conf1 == PTParams::numOfConf - 1)
			&& ((confsData_[conf1]).confDir == UP))
		(confsData_[conf1]).confDir = DOWN;
	if ((conf1 == 0) && ((confsData_[conf1]).confDir == DOWN))
		(confsData_[conf1]).confDir = UP;
	if ((conf2 == PTParams::numOfConf - 1)
			&& ((confsData_[conf2]).confDir == UP))
		(confsData_[conf2]).confDir = DOWN;
	if ((conf2 == 0) && ((confsData_[conf2]).confDir == DOWN))
		(confsData_[conf2]).confDir = UP;
}
void ParallelTempring::updateConfHistogram() {
	for (int i = 0; i < PTParams::numOfConf; i++) {
		if (confsData_[i].confDir == UP)
			histUp_[i]++;
		else
			histDown_[i]++;
	}
}
void ParallelTempring::initializeConfHistogram() {
	confData cd;
	std::fill(histDown_.begin(), histDown_.end(), 0);
	std::fill(histUp_.begin(), histUp_.end(), 0);
	std::fill(ptRate_.begin(), ptRate_.end(), 0);
	confsData_.clear();
	for (int i = 0; i < PTParams::numOfConf; i++) {
		cd.confNum = i;
		if (i % 2)
			cd.confDir = UP;
		else
			cd.confDir = DOWN;
		confsData_.push_back(cd);
	}
}
int ParallelTempring::getNumMeasure() {
	int min_num = QmcParams::minNumOfMeasure;
	for (int i = 0; i < PTParams::numOfConf; i++)
		min_num = fmin((confs_[i].getNumMeasure()), min_num);
	return min_num;
}
void ParallelTempring::printStatus() {
	static int precent = 0;
	double done = double(this->getNumMeasure())
			/ double(QmcParams::minNumOfMeasure);
	timeval currentTime, result;
	gettimeofday(&currentTime, NULL);
	timersub(&currentTime, &startAll_, &result);
	OutPutfilePT_ << "\n  Time :" << (double(result.tv_sec) / 3600.)
			<< " Num Of Closed : " << this->countClosed() << " Done :" << done
			<< "\n" << std::flush;
	this->printConf();
#ifdef SPLIT_SIM
	stopTimeEnd_ = (double(result.tv_sec) > PTParams::runTime);
	if (stopTimeEnd_)
	OutPutfilePT_ << "\n Time Ended" << endl;
#endif
}
void ParallelTempring::updateLambdas() {
	DoubleVec eta(PTParams::numOfConf - 1);
	DoubleVec delta_R(PTParams::numOfConf);
	DoubleVec f(PTParams::numOfConf);
	DoubleVec RSS = PTParams::RSS;
	double norm_eta = 0;
// calc delta_Lambda
	std::adjacent_difference(RSS.begin(), RSS.end(), delta_R.begin());
	delta_R.erase(delta_R.begin());
// Calc f
	std::transform(histUp_.begin(), histUp_.end(), histDown_.begin(), f.begin(),
			[](int u,int d)->double
			{	return double(u)/(double(d)+double(u));});
	std::adjacent_difference(f.begin(), f.end(), f.begin());
	f.erase(f.begin());
// Calc Eta
	for (int i = 0; i < PTParams::numOfConf - 1; i++) {
		eta[i] = sqrt(fabs(f[i])) / delta_R[i];
		norm_eta += sqrt(fabs(f[i]));
	}
	for (int i = 0; i < PTParams::numOfConf - 1; i++)
		eta[i] /= norm_eta;

// Calc New Rss
	for (int i = 1; i < PTParams::numOfConf - 1; i++) {
		double temp_sum = 0, target_sum = double(i)
				/ double(PTParams::numOfConf - 1), d_R;
		int j = 0;
		while (temp_sum < target_sum) {
			temp_sum += eta[j] * delta_R[j];
			j++;
		}
		j--;
		temp_sum -= eta[j] * delta_R[j];
		d_R = (target_sum - temp_sum) / eta[j];
		RSS[i] = PTParams::RSS[j] + d_R;
	}
	PTParams::RSS = RSS;
	Sigmas_.clear();
	Force_Factor.clear();
	for (int i = 0; i < PTParams::numOfConf; i++) {
		confs_[i].updateRS(PTParams::RSS[i]);
#ifdef ROT
		double lambda = 1
		/ (PI * double(QmcParams::NumP) * pow(PTParams::RSS[i], 2)
				/ (sqrt(3) / 2.));
#else
		double lambda = 1
				/ (PI * double(QmcParams::NumP) * pow(PTParams::RSS[i], 2));
#endif
		double epsilon = (1 / ((double) (QmcParams::M))) * QmcParams::beta;
		double sigma = sqrt(2 * lambda * epsilon);
		Sigmas_.push_back(0.5 * pow(sigma, -2));
		Force_Factor.push_back((1. / 9.) * pow(epsilon, 3) * (2. * lambda));
	}

	OutPutfilePT_ << "Lambdas: \n {" << std::flush;
	for (int i = 0; i < PTParams::numOfConf - 1; i++)
		OutPutfilePT_ << PTParams::RSS[i] << "," << std::flush;
	OutPutfilePT_ << PTParams::RSS[PTParams::numOfConf - 1] << "}\n"
			<< std::flush;
	this->maxEnergyDiff();
}
void ParallelTempring::adjustLamabdas() {
	long int run_num = 0;
	int runs = PTParams::numAdjustRuns;
	OutPutfilePT_ << "Start Adjust Lambda" << endl;
// First Run
//	for (int i = 0; i < PTParams::numOfConf; i++)
//		confs_[i].setNumOfThermal(QmcParams::numThermal * 10);
	this->thermalize();
//	for (int i = 0; i < PTParams::numOfConf; i++)
//		confs_[i].setNumOfThermal(QmcParams::numThermal);
///////////
	for (int j = 0; j < PTParams::numAdjustLambda; j++) {
		// Initialize
		int max_sweep_rate = this->getMaxRate() / 10;
		// run
		for (int run = 0; run < runs; run++) {
			this->runN(max_sweep_rate);
			for (int i = 1; i < PTParams::numAdjustPar; i++) {
				this->swapConfMove();
				this->updateConfHistogram();
			}
		}
		this->updateLambdas();
		this->dumpHistograms(j);
		this->initializeConfHistogram();
		runs *= 2;
		this->thermalize();
	}
	hsize_t dims = PTParams::numOfConf;
	hsize_t cdims = PTParams::numOfConf;
	String ds("DS1");
	String file_name = PTParams::outPathPT;
	file_name.append("RSS_Adjust.h5");
	hdf5_rw::dumpToHDF5Double(file_name, ds, &dims, 1, &cdims,
			&(PTParams::RSS[0]));
}
void ParallelTempring::dumpHistograms(int i) {
	hsize_t dims = PTParams::numOfConf;
	hsize_t cdims = PTParams::numOfConf;
	std::string ds("DS1");
	String file_name = PTParams::outPathPT;
	file_name.append("UpHist");
	if (i == -1)
		file_name.append("_Final");
	else
		file_name.append(boost::lexical_cast<std::string>(i));
	file_name.append(".h5");
	hdf5_rw::dumpToHDF5Int(file_name, ds, &dims, 1, &cdims, &(histUp_[0]));
	file_name = PTParams::outPathPT;
	file_name.append("DownHist");
	if (i == -1)
		file_name.append("_Final");
	else
		file_name.append(boost::lexical_cast<std::string>(i));
	file_name.append(".h5");
	hdf5_rw::dumpToHDF5Int(file_name, ds, &dims, 1, &cdims, &(histDown_[0]));
	file_name = PTParams::outPathPT;
	file_name.append("RSS");
	if (i == -1)
		file_name.append("_Final");
	else
		file_name.append(boost::lexical_cast<std::string>(i));
	file_name.append(".h5");
	hdf5_rw::dumpToHDF5Double(file_name, ds, &dims, 1, &cdims,
			&(PTParams::RSS[0]));
}
void ParallelTempring::maxEnergyDiff() {
	WLEnergyData E1, E2, E1Max, E2Max;
	double En1, En2, max_ef = 0;
	int i_max;
	for (int i = 0; i < PTParams::numOfConf; i++) {
		confs_[i].getEnergyCalcEnergy(E1, E2);
		En1 = E1.getEnergy();
		En2 = E2.getEnergy();
		if (fabs(En1 - En2) > max_ef) {
			max_ef = fabs(En1 - En2);
			E1Max = E1;
			E2Max = E2;
		}
	}
	OutPutfilePT_ << "Maximal Diff = " << max_ef << std::endl;
	OutPutfilePT_ << "Energy Calculated  = " << (E1Max.getEnergy())
			<< std::endl;
	OutPutfilePT_ << "Energy From WL Data  = " << (E2Max.getEnergy())
			<< std::endl;
}
//void ParallelTempring::print_graph(int precent) {
//	String FileNameSF(PTParams::Out_Path_PT);
//	FileNameSF += "SF_Graph";
//	FileNameSF += boost::lexical_cast<std::string>(precent);
//	FileNameSF += ".ps";
//	DoubleVec SF, WindSqr;
//	plstream *pls;
//	for (int i = 0; i < PTParams::numOfConf; i++) {
//		SF.push_back(confs_[i].get_SF_mean());
//		WindSqr.push_back(confs_[i].get_WindSqr_mean());
//	}
//	pls = new plstream();
//	// Initialize PLplot.
//	pls->star(1, 1);
//	// Select the multi-stroke font.
//	pls->fontld(1);
//	pls->sfnam(FileNameSF.c_str()); /* file name */
//	//	pls->sdev("png"); /* device type */
//	pls->init();
//	pls->env(Lambdas_.back(), Lambdas_.front(), 0, 0.3, 0, 0);
//	pls->lab("Lambda", "SF", "SF vs Lambda");
//	pls->line(Lambdas_.size(), &(Lambdas_.at(0)), &(SF.at(0)));
//	delete pls;

